Theory of vortex-lattice melting in a one-dimensional optical lattice 
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We investigate quantum and temperature fluctuations of a vortex lattice in a one-dimensional 
optical lattice. We discuss in particular the Bloch bands of the Tkachenko modes and calculate the 
correlation function of the vortex positions along the direction of the optical lattice. Because of 
the small number of particles in the pancake Bose-Einstein condensates at every site of the optical 
lattice, finite-size effects become very important. Moreover, the fluctuations in the vortex positions 
are inhomogeneous due to the inhomogeneous density. As a result, the melting of the lattice occurs 
from the outside inwards. However, tunneling between neighboring pancakes substantially reduces 
the inhomogeneity as well as the size of the fluctuations. On the other hand, nonzero temperatures 
increase the size of the fluctuations dramatically. We calculate the crossover temperature from 
quantum melting to classical melting. We also investigate melting in the presence of a quartic 
radial potential, where a liquid can form in the center instead of at the outer edge of the pancake 
Bose-Einstein condensates. 

PACS numbers: 03.75.Lm, 32.80.Pj, 67.40.-w, 67.40. Vs 



I. INTRODUCTION 

Since the onset of experiments on Bose-Einstein con- 
densates, vortices have attracted a lot of attention. When 
the Bose-Einstein condensate is rotated faster than some 
critical rotation frequency Q c , a vortex appears in the 
gas Ql Hj. Q Upon increasing the rotation frequency fur- 
ther, the number of vortices increases 0, IE IS j ancl 
the vortices order themselves in a hexagonal Abrikosov 
lattice H [H [H d d Q If the rotation fre- 

quency is increased even further, the very rapidly rotat- 
ing ultracold bosonic gases have been predicted to form 
highly-correlated quantum states [H [H 13 EES In 
these states, the Bose-Einstein condensate has been com- 
pletely depleted by quantum fluctuations, and quantum 
liquids appear with excitations that can carry fractional 
statistics. Some of these states have been identified with 
(bosonic) fractional quantum Hall states [T^.l2lll2^l. such 
as the Laughlin state [23J , the Moore- Read state |24j and 
various Read-Rezayi states pM l2r| . 

In this article we study the quantum fluctuations of 
vortices in a one-dimensional optical lattice. Optical lat- 
tices provide a powerful tool in manipulating a Bose- 
Einstein condensate. By using a three-dimensional op- 
tical lattice, the Bose-Hubbard model j2]j was experi- 
mentally realized and the superfluid-Mott insulator tran- 
sition was observed [2^. Moreover, a two-dimensional 
optical lattice was used to create one-dimensional Bosc- 
gases and to study the crossover between a superfluid 
Bose gas and the "fermionized" Tonks gas |2|J, |3(j . The- 
oretically, the combination of a two-dimensional optical 
lattice and a single vortex was predicted to give inter- 
esting effects around the superfluid-Mott insulator tran- 
sition [3l|]. Optical lattices can also be used to provide 
a pinning potential for vortices [32l I33I l34l I3M l3q , and 
experiments on this topic are ongoing |37j |. The physics 
of a single vortex line in a one-dimensional optical lattice 
is recently extensively studied [H |H EE E3 IH E3 . By 



putting fcrmions in the vortex core, this system can even 
be used to create a superstring in the laboratory 01 ■ 

A one-dimensional optical lattice divides the Bose- 
Einstein condensate into a stack of two-dimensional pan- 
cake condensates that are weakly coupled by tunneling 
as schematically shown in Fig. ^ An important con- 
sequence of this setup is that the modes of the on-site 
two-dimensional vortex lattice form Bloch bands as a 
function of the axial momentum. In this paper we pay 
special attention to the Tkachenko modes [4f| , which are 
the lowest-lying modes of the vortex lattice. Recently, 
these modes have been investigated for a single two- 
dimensional vortex lattice both exp erimentally [T3 . IT^j 
and theoretically [MEllSIMIilllJIIllII. 

The number of particles in a single pancake is much 
smaller than in a Bose-Einstein condensate in a harmonic 
trap, and hence the fluctuations are much larger. There- 
fore, this system is a promising candidate to reach the 
quantum Hall regime. The requirement for that is that 
the ratio v = N/N v of the number of atoms N and the 
number of vortices N v is smaller than a critical value v c . 
The ratio v plays the role of the filling factor and esti- 
mates for the critical v c are typically around 8 [T^. l54j . 
However, observed filling factors arc up till now always 
greater than 100, where almost perfect hexagonal lattices 
form and no sign of melting can be seen [lj|. These ex- 
periments are carried out with Bose-Einstein condensates 
consisting of typically 10 5 particles, whereas the maxi- 
mum number of vortices observed is around 300. De- 
creasing the number of particles results in loss of signal, 
whereas the number of vortices is limited by the rotation 
frequency that has to be smaller than the transverse trap- 
ping frequency. Adding a quartic potential, which stabi- 
lizes the condensate also for rotation frequencies higher 
than the transverse tra p fr equency, has until now not 
improved this situation |l5| . although it has opened up 
the possibility of forming a giant vortex in the center 
of the cloud [KM l56l l57j . Applying a one-dimensional 




II. EFFECTIVE LAGRANGIAN 

In this section we derive the effective lagrangian for the 
vortex fluctuations. Using the ansatz that the conden- 
sates wavcfunctions are part of the lowest Landau level, 
we find the classical groundstates and expand the action 
up to second order in the fluctuations. 



A. Tight-binding approximation 



FIG. 1: (Color online) Setup in which the melting of the 
vortex lattice is studied. An optical lattice along the z direc- 
tion divides the condensate into pancake condensate which 
are coupled by tunneling processes. 



optical lattice produces pancake condensates each con- 
taining typically 1000 particles, such that quantum fluc- 
tuations are strongly enhanced. Experimental signal is 
still conserved because of the combined signal of all the 
pancake. Moreover, because of the small number of par- 
ticles in each pancake shaped Bosc-Einstcin condensate, 
finite-size effects become very pronounced in this setup. 
In particular, the critical filling factor for the melting 
of the lattice v c changes compared to the homogeneous 
situation. As a further consequence, melting is not ex- 
pected to occur homogeneously but starts at the outside 
and then gradually moves inwards as the rotation speed 
increases |5||- Therefore, phase coexistence is expected, 
where a vortex crystal is surrounded by a vortex liquid. 

The optical lattice gives also the exciting possibil- 
ity to study the dimensional crossover between two- 
dimensional melting and three-dimensional melting, by 
varying the coupling between the pancakes. This system 
also exhibits an intruiging similarity to the layered struc- 
ture of the high-T c superconductors |5^, |6(| • Recently, 
the density profiles for quantum Hall liquids in this ge- 
ometry have been calculated Also the static prop- 
erties of the lattice in a double and multilayer geometry 
have been investigated E3 • Finally, classical melt- 
ing between shells of vortices in a single two-dimensional 
Bose-Einstein condensate has also been studied recently 

In this article we study this interesting physics by ex- 
panding on our previous work j6^. The paper is or- 
ganized as follows. In Sec. El we present the theoreti- 
cal foundations of our work and derive the effective la- 
grangian for the vortex fluctuations. In Sec. IIIII we dis- 
cuss the excitation spectrum and we pay in particular 
attention to the Bloch bands of the Tkachenko modes. 
In Sec. IIVI we present result on the quantum melting, 
and in Sec. [V] we take into account the effects of a 
nonzero temperature. We calculate the crossover tem- 
perature from quantum melting to classical melting. In 
Sec. I VII we present results on the addition of a quartic 
potential. We end up with conclusions in Sec. IVIII 



We closely follow the approach in Ref. |4jj for a single 
vortex, extending it to the case of a vortex lattice. Our 
starting point is a cigar-shaped Bosc-Einstcin condensate 
trapped by the potential 



nr) = ^>ir 2 +^z 2 ), 



(1) 



where m is the atomic mass and oj r and oj z are the radial 
and axial trapping frequencies, respectively. Since we 
assume a cigar-shaped trap, we have that u> z <C uj±_ and 
we neglect the axial trapping in the rest of this work. 
Experimentally this can also be realized by using end- 
cap lasers, which results in a box-like trapping potential 
in axial direction (6(|. In addition, the Bose-Einstein 
condensate experiences a one-dimensional optical lattice 



V z (r) =V 2 sin 2 



2ttz 



(2) 



where V z is the lattice depth and A is the wavelength of 
the laser light. The lattice potential splits the conden- 
sate into N s two-dimensional condensates with a pan- 
cake shape. Each of the pancake Bose-Einstein conden- 
sates contains on average N atoms, so the total number 
of atoms is N S N. We consider a deep optical lattice, 
such that its depth is larger than the chemical poten- 
tial of the two-dimensional condensate, but still we as- 
sume that there is full coherence across the condensate 
array. This means in particular that the lattice potential 
should not be so deep as to induce a superfluid-Mott in- 
sulator transition. Typically the required lattice depth to 
reach the supcrfluid-Mott insulator transition in a three- 
dimensional lattice with one particle per lattice site is of 
the order of 10-E 2 , where E z is the recoil energy of an 
atom after absorbing a photon from the laser beam and 
given by 



E, 



(2irk/\) 2 
2m 



(3) 



In a onc-dimcnsional lattice the number of atoms in 
each lattice site is typically much larger than in a three- 
dimensional lattice and the transition into the insulating 
state requires a much deeper lattice |(37l ■ The superfluid- 
Mott insulator transition in a one-dimensional optical 
lattice has recently been observed [6^, but also in that 
case the number of particles per lattice site is of the order 
of one. 
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The action describing the system in the rotating frame 
is given by S = J dt J d 3 x £(x, t), where the Lagrange 
density is given by 

m + - Kx(x) + ql z - ||*| 2 ) *. 

(4) 

Here \I/(x, t) is the Bosc-Einstcin condensate wavefunc- 
tion, and f2 is the rotation frequency. Moreover, 



L- 



ih(yd x - xdy) 



(5) 



is the angular momentum operator, and g = 4irh 2 a/m 
is the interaction strength, with a the three-dimensional 
scattering length. 

Since we assume a deep optical lattice potential, we 
can perform a tight-binding approximation and write 



(6) 



where i labels the lattice sites and z% = iX/2 is the posi- 
tion of the ith site. The wavefunction in the z direction 
&(z) is chosen to be the lowest Wannier function of the 
optical lattice. Because of the deep optical lattice this 
wavefunction is well approximated by the ground-state 
wavefunction of the harmonic approximation to the lat- 
tice potential near the lattice minimum. The frequency 
associated with this harmonic trap is 



2?r 2V z 



A V m 

and the wavefunction <J>(z) is given by 



$(z) = (nil) 1 / 4 



exp 



2t*J' 



(7) 



(8) 



where II — \J hjrauj^. 

Writing the Lagrange density as C = T — £, the time- 
derivative part of the Langrage density becomes 



iH 



(9) 



Terms coupling wavefunctions on neighboring sites do not 
appear in this kinetic part of the lagrangian, because of 
the orthogonality of the Wannier functions on different 
sites. The energy part of the Lagrange density can be 
written as 



(10) 



with 

Si = 



2m 2 



Jij = -t + , 



(11) 



and (ij) denotes that the sum is taken over nearest- 
neighboring sites. We have defined the effective two- 
dimensional coupling strength 



g' = g J dzMztf 



Airfr 2 



2ir£j,m 



(12) 



and the hopping amplitude 

h 2 d 2 



t 



d 2 *(z) 



2m 



+ V z {z) 



*(z + A/2). (13) 



Using the Gaussian wavefunction from Eq. © to calcu- 
late this hopping amplitude underestimates the hopping 
amplitude, since this approximation is only good in the 
vicinity of the center of the wells of the optical lattice, 
and not in the classically forbidden regions where the 
overlap between the neighboring wavefunctions is maxi- 
mal. Therefore, we use the result 



t = AV^E^ 4 exp[- 




(14) 



which is exact for a deep optical lattice |69| . The parame- 
ters g' and t are both fully determined by the microscopic 
details of the atoms and the optical lattice. 



B. Lowest Landau level approximation 

Melting is only expected for a Bose-Einstein conden- 
sate that is weakly interacting. The transverse wavefunc- 
tion can then be taken to be part of the lowest Landau 
level |73|. This implies that we describe the Bose- 
Eintcin condensate as a compressible fluid. Thus we con- 
sider wavefunctions of the form 

y, t) « ]J{w - w m (t)) exp[-M72], (15) 



where iu = (x + iy)/£, w ni (t) = (x ni (t) + iy ni (t))/£ and 
x m(*) = {x n i{t),yni(t)) is the position of the n th vortex 
at site i. The vortex positions are the dynamical vari- 
ables and are therefore time dependent. Here £ is the 
"magnetic length" , which is normally identified with the 
radial harmonic length £± = \J h/mui±. The validity of 
the lowest Landau level approximation can then in mean- 
field theory be estimated to be limited by the condition 
III 



VtoN— < ( 1 - — ) , 



(16) 



which is recently compared with exact diagonalization 
results for small numbers of particles |7^| . To increase the 
validity of this study, we use £ as a variational parameter 
instead of fixed it to the radial harmonic length, such that 
our results are also valid for stronger interactions and 
slower rotation [j4|. The associated frequency is 10 = 



4 



h/m£ 2 . In practice, it turns out that to is always well 
approximated by the rotation frequency CI. 

Making use of the fact that the wavefunctions <j> n (w) = 
w n e -\w\ ji j w £2 n \ f orm a complete and orthonormal ba- 
sis for the lowest Landau level wavefunctions, it is easy 
to derive that within these approximations we have that 



d 2 r *t 



2m 



® i= 2~P / d2rr2n '( r )' 



(17) 



where the density rij(r) is a function of the vortex posi- 
tions X„i(t). 

From now on distances are rescalcd by £, frequencies 
are scaled by the radial trapping frequency u>j_, and we 
define a dimensionless interaction strength by means of 



U = N H^ = N . 



2ir£, 



(18) 



The on-site part of the energy functional can then be 
written as: 

j^jf = \ (w + ~ - 2fl) r 2 n t (r) + 2™Un 2 {v). (19) 

The lowest Landau level wavefunctions and, therefore, 
also the atomic density, are fully determined by the lo- 
cation of the vortices. To consider the quantum mechan- 
ics of the vortex lattice we, therefore, replace the func- 
tional integral over the condensate wavefunctions by a 
path integral over the vortex positions. This involves 
a non-trivial Jacobian, which does not change the re- 
sults presented here, because we always consider the case 
that N 3> 1. In the calculations we take the scattering 
length of 87 Rb, A = 700 nm and V z /E z = 16, which 
gives U = 25N. The qualitative features, however, do 
not depend on the value of U. 



C. Classical Abrikosov lattice 

To determine the quadratic fluctuations around the 
Abrikosov lattice, we first have to find the classical 
groundstate. We calculated this groundstate for up to 
37 vortices. The number of vortices in the condensate 
increases with the rotation frequency. For small num- 
bers of vortices, the groundstate is distorted from the 
hexagonal lattice [TBI [7(| . In general, there are also vor- 
tices far outside the condensate. When there arc more 
than 18 vortices in the condensate, there is generally one 
vortex in the center, while the other vortices order them- 
selves in rings of multiples of six vortices. However, when 
a new ring of vortices enters the condensate, there is an 
instability towards an elliptic shape deformation. This 
is related to the elliptic shape deformation that occurs 
before a single vortex enters the condensate , an d that 
has been investigated before theoretically (z3> EE EH 
and has also been observed [8l]]. This shape instabil- 
ity plays an important role in the classical melting of 



the vortex lattice, as we show lateron. Pictures of the 
classical groundstate are given in Fig. [21 for fixed inter- 
action U and different rotation frequencies. The number 
of vortices within the condensate as a function of rota- 
tion frequency is plotted in Fig. where we limited our 
study to vortex lattices that exhibit hexagonal symme- 
try around the origin. In the fluctuation calculation wc 
only consider the regime where configurations consisting 
of one vortex in the center surrounded by rings of six vor- 
tices are stable. In that case, the coarse-grained atomic 
density is well approximated by a Thomas-Fermi profile 



D. Fluctuations 

Next we study the quadratic fluctuations by expanding 
the action up to second order in the fluctuations u„j = 
x„i — (x n i). This yields an action of the form 



S = ■ {Tid t - E) ■ Ui - tJ2 u * ' J ' u j 



(20) 



where u» = (. . . , u m ;, . . .) is the total displacement vector 
of all the point vortices on site i. The matrices T, E 
and J depend on CI, U, and the classical lattice positions 
(x„j), and are formally given by 



E,, 

T„ 
J,, 



E x 

En 



= d 2 v 



= d 2 Y 



d 2 r 
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dXni 





Jit 



(21) 



We also expand the density in the fluctuations by means 
of 



ni(r) 



( r ) ■ u ii 



(22) 

where no (r) is the equilibrium particle density associated 
with the classical lattice. Wc have defined the vectors 



n„(r) 



n x„ (r) 



(23) 



which form the dipolc densities and are associated with 
linear variations around the classical solution, and the 
tensors 



n nm (r) 



n Xn x m (r) n Xn y m (r) 
n yn x m (r) n ynVm (r) J ' 



(24) 



which form the quadrupole densities and are associated 
with the quadratic variations around the classical solu- 
tion. Moreover, 



n„ = J d 2 r n„(r) 



(25) 
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FIG. 2: (Color online) Classical vortex lattice and density profile for rotation frequencies Q/lj± — .8, .91, .93 .95 and .97. Here 
U — 10, which corresponds to N — 250. White means high density, black low density. The vortex positions are indicated by a 
dot, such that also the vortices outside the condensate are visible. 
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FIG. 3: The number of vortices N v that is within the con- 
densate as a function of the rotation frequency Q. Depending 
on the rotation frequency, the groundstate is either a vor- 
tex lattice with a vortex in the center surrounded by rings 
of a multiple of six vortices such that the number of vortices 
can be written as 1 + 6n r , or a configuration where there is 
no vortex in the center and the vortices order themselves in 
rings of multiples of three vortices, such that the total num- 
ber of vortices in the condensate is a multiple of three. For 
Q/c^x > -9, the first situation is the groundstate (except for 
a small region, where the Bose-Einstein condensate becomes 
elliptically deformed) and we only calculate the fluctuations 
for this case. For this plot the condensate radius was defined 
as the radius where the angular averaged density drops below 
0.003. 



and 



d 2 r n nm (r), 



(26) 



such that Tij(r) is always normalized to 1. The density 
no(r), and the tensors n„(r), n nm (r) are independent of 



time and the layer index i, since only the fluctuations 
in the vortex positions Uj(i) are taken as the dynamic 
variables. Making use of the fact that 



th(t) oc Y[\w - Wj| 2 exp[-|K;| 2 ], 



(27) 



the following expressions can be derived for the dipole 
density 



/ n 1\X \X n i)) . . 

= |r-(x M )P "°< r >» 

/ x -% - (ym)) , , 

n »»( r ) = \7ZT7Z AI2 n o( r )- 



(28) 



For the quadrupole density for a single vortex we get 
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whereas for two different vortices we find 

n>x nVm (r) 
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Mr), (30) 
n (r), 
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All matrices in the action for the fluctuations can be 
completely expressed in terms of these functions. From 
Eq. lfT9")l we read off that 
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E m „ = / <rr \ -{u + 1/cj - 2ft)r 2 |^(n„n m - n nm )n (r) - -(n„n m (r) + n m n„(r)) + n m „(r) 



(31) 



+2itujU 



(3n„n m - 2n„ m rio(r) - 2 (n„n m (r)n (r) + n m n„(r)n (r)) + n„(r)n m (r) + 2n„ TO (r)n (r) 



From a straightforward derivation we obtain also 
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(33) 



The calculation is simplified considerably by making use 
of the hexagonal symmetry of the vortex lattice. Because 
of this symmetry the number of matrix elements to be 
calculated for each matrix is less then 12N% + 32N V in- 
stead of 72N%, when all matrix elements are calculated 
independently. 
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FIG. 4: Dispersion of the Tkachenko modes along the direc- 
tion of the optical lattice. For this plot the parameters were 
chosen as: U = 10, t = 1/10, and Q — .97. One gapless linear 
mode and five gapped tight-binding-like modes can be iden- 
tified. The gapless mode is the acoustic Goldstone mode as- 
sociated with the broken 0(2)-symmetry due to the presence 
of the vortex lattice. Note that these results were obtained in 
the lowest Landau level approximation, which corresponds to 
the compressible limit. 



E. Diagonalization 

To diagonalizc this action along the z axis, we perform 
a Fourier transformation to obtain 

S = u *k ■ i Tid t — E — t[l — cos(fcA/2)J]) ■ u fc . (34) 



Finally, we completely diagonalize this action by a trans- 
formation 

v fe = P k u k , (35) 
that is normalized such that the action becomes 

S = ^2,vl a (id t - oj a (k))v ka , (36) 

k . a 

where to a (k) are the mode frequencies of the vortex lat- 
tice. This means that the v ka , where k labels the mo- 
mentum in the z direction and a labels the mode, cor- 
respond to bosonic operators with commutation relation 
[vka, Vk'a' \ = Skk'Saa'- This allows us to calculate the 
expectation value for the fluctuations in the vortex po- 
sitions, but also for the correlations between the various 
point vortices. 



III. TKACHENKO MODES 

The Tkachenko modes are almost transverse modes of 
the vortex lattice. In a harmonic trap with cylindrical 
symmetry they become modes which are almost angu- 
lar. In the radial direction their spectrum is discretized, 
because of the finite lattice size. The number of radial 
Tkachenko modes equals the number of rings in which 
the vortices have ordered themselves. For 37 vortices 6 
Tkachenko modes can be identified. A close comparison 
with continuum theory for a finite-size system, where also 
a discrete spectrum was found , is possible but beyond 
the scope of this work. Moreover, the Tkachenko modes 
also have a dispersion in the z direction. Without the op- 
tical lattice some aspects of these modes were recently in- 
vestigated |83| . For typical parameters this dispersion is 
plotted in Fig. 0] while the modes are displayed in Fig. 
As can be seen in the latter figure, some of the vortices at 
the edge deviate from having a purely angular motion. As 
is clearly visible, there is one gapless mode, which is linear 
at long wavelengths, while the other modes are roughly 
just tight-binding-like. Moreover, various avoided cross- 
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FIG. 5: (Color online) Tkachenko modes. For this plot the 
parameters were chosen as: U = 10, O = .97, and k — 0. 
Mode 1 is the Goldstone mode and corresponds to a pure 
rotation, while the modes 2-6 are gapped and ordered along 
increasi ng g ap. They were previously called s-band modes 
in Ref. [l3jl . As is visible, some of the vortices at the edge 
deviate from having a purely angular motion. 




FIG. 6: Node structure of the Tkachenko modes. The modes 
are ordered in the same way as in Fig. |S] Points are connected 
by straight lines to increase the visibility. When a wavelength 
is associated with the position of the nodes, it can be seen 
that for increasing gap the wavelength becomes shorter. This 
agrees qualitatively with the continuum theory of Ref. IH^l . 



ings between these modes are clearly visible. The gaplcss 
mode is the Goldstone mode associated with the spon- 
taneously broken rotational 0(2)-symmctry due to the 
presence of the vortex lattice. When the tunneling rate 
is very small, the gapped modes have exactly a tight- 
binding dispersion and the gapless mode gets a disper- 
sion proportional to sin(fcA/4). This can be understood 
by observing that in this case the modes are decoupled 
and the hamiltonian for the gapless mode reduces to the 
Joscphson hamiltonian 



H 



Ej cos( 

(ij) 



(37) 



Writing = —idtf,, the action in momentum space thus 
reads after a quadratic expansion of the Josephson energy 



S = J 



2[1 - cos(k\/2)]Ej -ioj 
iuj Ec 



from which we deduce the dispersion 



(38) 



oj 2 (k) = 2[1 - cos(k\/2)}EjE c = 4:Sm 2 (k\/4)EjE C - 

(39) 

It is interesting to note that for a small rotation 
frequency, which implies a small vortex lattice, the 
Tkachenko modes are not the lowest-lying modes. For 
U = 10, a Tkachenko mode becomes the lowest-lying 
gapped mode when fl > 0.978, but there are many modes 
in between the second and the third Tkachenko mode. 
This confirms the expectation that increasing the vortex 
lattice will bring down the Tkachenko spectrum more and 
more. 



IV. QUANTUM MELTING 

In this section we study vortex-lattice melting due to 
quantum fluctuations. We apply the Lindcmann crite- 
rion to estimate the position of the melting transition. 
We study the influence of tunneling between the pancake 
condensates and compare with a local-density approxi- 
mation. By looking at correlations between the vortices, 
we can distinguish between various phases, which are par- 
tially melted. 



A. Single-layer geometry 

Quantum fluctuations of the vortices ultimately result 
in melting of the vortex lattice. To decide whether or 
not the lattice is melted, we use the Lindemann criterion, 
which in this inhomogeneous situation has to be applied 
locally. The Lindemann criterion means that the lattice 
is melted, when 
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FIG. 7: Crystal radius R CI normalized to the condensate ra- 
dius R as a function of the rotation frequency for N = 250 
and U = 10. The solid line is for t = and the dashed lines 
are for t = 1/1000 and 1/100 (from left to right). The dotted 
lines are the result of the local-density approximation, where 
for the upper curve we took the criterion N/N v — 6 and for 
the lower curve N/N v — 8. In all the figures in this paper we 
have connected the points by straight lines. Due to the finite 
numerical resolution the jump between the plateau's appears 
therefore not vertical. 



where A n i is average distance to the neighboring vortices. 
The critical value cl, is known as the Lindcmann pa- 
rameter. This parameter is a priori unknown, but values 
found from comparison with Monte Carlo simulations are 
typically in the range cl = 0.1 — 0.3. The value cl = 0.1 
was recently found from elaborate calculations for a tri- 
angular vortex lattice in high-T c superconductors, which 
also compared well to experiments in that case [84|. We, 
therefore, use this value in our calculations. Note that 
the results we obtain depend quantitatively on the value 
of the Lindemann parameter. Changing this value shifts 
the curves, but the qualitative features remain the same. 

Due to the inhomogeneity of our system, we have to 
apply this criterion locally. Because the coarse-grained 
particle density decreases with the distance to the origin, 
vortices on the outside are already melted, while the in- 
ner part of the crystal remains solid j5^- Therefore, a 
crystal phase in the inside coexists with a liquid phase 
on the outside. In Fig. 0we compute the radius of the 
crystal phase R CI normalized to the condensate radius R, 
as a function of the rotation frequency for fixed a num- 
ber of particles and a fixed interaction strength U, but 
for various hopping strengths t. Also here we define the 
condensate radius R as the radius for which the angu- 
larly averaged density drops below 0.003. The crystal 
radius R cr is defined as the radius of the innermost ring 
of vortices that is melted according to the Lindemann 
criterion. When according to this definition R cr > R, 
we set the crystal radius equal to the condensate radius, 
i.e., R CI = R. The ratio R C r/R shows discrete steps be- 
cause of the ring-like structure in which the vortices order 



B. Local-density approximation 

We compare this with a simple local density calcula- 
tion, where the criterion N/N„ = n(r)/n v (r) = 6 |l9l | or 
N/N v = n{r)/n v {r) = 8 [54J is applied locally, by mak- 
ing use of a Thomas-Fermi density profile to describe the 
coarse-grained atomic density. Substituting the Thomas- 
Fermi profile 



n T F(r) 



ttR\ f 
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p2 

n, TF 



(41) 



where -Rtf is the Thomas-Fermi radius, in the on-site 
energy in Eq. (|19|) . we obtain 
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We have introduced here also the Abrikosov parameter 
(3 = j n 2 (r)/(J n(r)) 2 ~ 1.1596 for the hexagonal lattice. 
Minimizing for Rtf gives 



16(3Ulu 2 



TF 



i + to 2 - 2ujn ' 

which on substitution gives for the on-site energy 



(42) 



This is minimized for 



0, 



which sets the variational parameter 
Fermi radius is then given by 



(43) 

The Thomas- 
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TF 
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n 2 /3u 



i-n 2 ' 

In lowest order the vortex density is given by 

1 



n v r 



Solving now 



j/(r) 



n(r) 
n v (r) 



(44) 



(45) 



(46) 



we get for the crystal radius normalized to the Thomas 
Fermi radius the result 



-Rcr 

Rtf 



v c R 2 
2N 



\ 




(47) 



For comparison this line is plotted in Fig. [7| As can 
be seen there are important finite-size corrections to the 
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local-density approximation. The local-density approxi- 
mation fails to take into account the discrete nature of the 
vortex positions. Moreover, it predicts the total melting 
of the crystal at considerable higher rotation frequencies 
than the exact answer. 

In principle there are corrections to the vortex density 
due to the Thomas-Fermi profile of the particle density. 
A better approximation near the center of the trap is 



>ity. 
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(48) 



However, this vortex density becomes zero well within the 
Thomas- Fermi radius. As a result the ratio n(r)/n v (r) 
diverges and is always bigger than 8. Higher-order con- 
tributes should be taken into account to solve this prob- 
lem. Instead we compare with the constant vortex den- 
sity 



(49) 



TF 



which corresponds to a Gaussian distribution of the 
atomic density with radius i?TF- Straightforward deriva- 
tion gives then 



Rtf 
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2N 
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n 2 pu 



(50) 



(51) 



This shifts the curves for the local-density theory in Fig. 
a little bit upwards, and makes the comparison even 
less favourable. 
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FIG. 8: Crystal radius R CI normalized to the condensate ra- 
dius R as a function of the hopping amplitude t for TV = 250, 



U = 10 and ST2 



.96. 
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C. Multi-layer geometry 

When the tunneling between pancakes is turned on, the 
fluctuations are also coupled in the axial direction. This 
decreases the fluctuations in the vortex displacements be- 
cause the stiffness of the vortices increases. Hence melt- 
ing occurs for higher rotation frequencies, as is visible in 
Fig. [7| In Fig. [S]we show the crystal radius for fixed ro- 
tation frequency and increasing hopping amplitude. To 
determine the presence of crystalline order in the axial 
direction, we calculate the correlation function 

^ e tq.(Uni-Un i )) ) 

which is related to the static structure factor. For our 
gaussian theory, this reduces to 

e -<[q.(u„,-u„ J )] 2 )/2^ 

For the central vortex the correlation function (uoiUfy) 
decays exponentially, which signals long-range order. For 



FIG. 9: Axial correlation function e -«<i-( u «i- u »i)) >/ 2 as a 
function of the distance along the axial direction for q = 
(1/1,1/1), N = 250, U = 10, and fi = 0.97. The vari- 
ous curves have, from top to bottom, a hopping parameter 
t = 10~ 2 , t = 10~ 3 , t = 10 -4 , t = 10 -5 , and t = 10 -6 , re- 
spectively. The correlation function decays as a power law 
and in the inset the power is plotted as a function of the 
hopping amplitude t. The power thus scales with \/\ft. 

the other vortices ([q- (u n j — u„j)] 2 ) grows as a logarithm, 
such that e~^ q '( Urei-u " J )l 2 ^ decays algebraically. This is 
in agreement with the expectation for a one-dimensional 
system at zero temperature. In Fig. Elwc show the axial 
correlation function for one of the vortices outside the 
center. 

Since the lattice spacing A/2 is the only length scale 
in the z direction, it also determines the axial correlation 
length. The power of the algebraic decay scales with 
1 / y/i. This can be understood when we assume that the 
gapless Tkachcnko mode gives the dominant contribution 
to the axial correlation. Using the effective hamiltonian 
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for this mode from Eq. I|37|) . wc obtain the following 
expression: 



2 



—==±== x 

v EjEc 

Ec 



dk sin 2 (fc(i - j)/2)) 



sin(fc/2) 



(52) 

Ejsm(k/2)] , 



where we rescaled k , in units of the lattice spacing A/2, 
r„ is the distance of the vortex to the origin and we used 
that 1 - cos(fc(i - j)) = 2sin 2 (fc(i - j)/2). Only the 
first term in the brackets is divergent for k — > 0, so wc 
neglect the second term. We make the approximation 
sin(fc/2) ~ k/2. When i — j is large we can approximate 
the rapidly oscillating function sm 2 (k(i — j)/2) by its 
average 1/2, except for a region near k = 0. We obtain 
then 



sm 2 (k(i-j)/2) 



k 

sin 2 (fc(z-j)/2) 




(53) 



2 J JTT. k 

dk— + -Iog4(»- j) 

2r 2 /Er 

-aj-g log4 i-i +d . 

7T V £j 

where the constant Ci = 7 + log(87r) + Ci(87r) ~ 3.80296. 
Using now the fact that Ej oc t, wc conclude that indeed 
the power of the algebraic decay of e _ ^ q ( - Uni_Unj -" ^ 2 
scales with 1 / \/t. 

For nonzero temperatures we obtain a linear rise of the 
correlation function ((u„i — u n j) 2 ), which can be under- 
stood from the same argument, since then we have to 
add the Bosc-Einstcin distribution, which for very low 
temperatures can be approximated by 



1 



k H T 



k R T 



J (k)/k B T_ 1 2^E.jE c sm(kX/4)' 



(54) 



When we again rescale the momentum by the lattice 
constant A/2 and approximate sinfc/2 ~ k/2, we ob- 
serve that the 1/fc-factor in the intergrand of Eq. I|53|) 
is replaced by l/k 2 . This gives that up to a constant 
((uj„ - Uj n ) 2 ) oc i — j and e~^ q '( Ur "~ u " J )) ^ 2 decays ex- 
ponentially. 



D. Correlations 

Melting for small arrays of electrons in quantum dots 
that are composed of two or three shells [8J, |83, |83, |88j , 
and of vortex lattice shells 0, is predicted to occur 
in two stages. First the oriental order between different 
shells is destroyed, while after that the radial order is de- 
stroyed. We also find this behavior for the vortex lattice. 
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FIG. 10: (Color online) Phases that can be distinguished by 
looking to radial and angular correlations between vortices on 
the same and on different shells: a) a solid, b) a solid of solid 
shells, c) a hexatic liquid, d) a liquid. 



For that we compute the correlation function between 
various vortices and decompose it in radial and angular 
fluctuations. A shell structure is found, where vortices 
with almost (but not necessarily exactly) the same dis- 
tance group together. Between these shells angular fluc- 
tuations dominate, while within the shell angular fluctu- 
ations are suppressed and radial fluctuations dominate. 
This can be understood if we assume that the Tkachenko 
modes dominate the fluctuations, because they leave the 
rings intact, but change the relative angle between the 
rings. To quantify this, wc introduce the following order 
parameter to measure fluctuations in the relative angle 
of neighboring shells 



(55) 



for neighboring vortices n and m on different shells, 
where r n is the distance of vortex n to the origin and 
only the angular part of the fluctuations in the displace- 
ment field is considered. Let 7V™ ax be the maximum 
number of vortices on either of the two shells. The shells 
arc decoupled with respect to each other if 



r SS > 

nm 



2tt 

AT ma: 
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(56) 



with cl = 0.1 the same Lindemann parameter. Further- 
more we introduce the correlation function 



r — 

- 1 - n.m. — 



2A 2 

rnn 



(57) 



where A m „ is the distance between the neighboring vor- 
tices m and n and we split this up in a radial and an- 
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FIG. 11: Phase boundaries between the solid (S), the solid 
of solid shells (SS), and the hexatic liquid (HL) as a function 
of radial distance and rotation frequency for N = 250 and 
U = 10. The solid line is for t = and the dashed line is for 
t = 1/100. 



gular component. Moreover we distinguish between the 
case when n and m are on the same shell, which we de- 
note by T T r ad and r^ ng , or on different shells, denoted by 
T^f and r™f. We use the definition that order is de- 
stroyed if the order parameters exceeds the Lindemann 
parameter, i.e., r > c| = 0.01. When the fluctuations 
are small, the vortex lattice is in the crystal phase. This 
phase is defined as T nm = + T^f < c| , for n and m 
on the same shell and on different shells. In particular, 
there is positional order of the radii of the different shells. 
However, we can distinguish between the case that there 
exists oriental order between the shells, in which case 
the relative angle is locked and r^ n < c 2 L . This 
we call the solid phase (S). When the oriental order be- 
tween the shells is destroyed we call this phase the solid 
of solid shells (SS). The melted phase is charactarized 
by T nm = T 1 ^ + 1^"! > c| , and in particular the posi- 
tional order in the radii of the different shells is destroyed. 
Within this phase it is still possible to have a well-defined 
angle between vortices on the same shell, i.e., r^ ng < c\. 
This we call the hexatic liquid (HL), whereas when this 
order is destroyed we call it the liquid (L). These phases 
are schematically indicated in Fig. ^| 

These criteria should again be applied locally because 
of the inhomogeneous density, and we define the phase 
boundaries in analogy to the previous definition as the ra- 
dius of the innermost ring that is partly melted. For our 
parameters it turns out that the angle between neighbor- 
ing vortices on the same ring is always well-defined, but 
we can identify the transition between the solid, the solid 
of solid shells, and the hexatic liquid. The result of this 
calculation is presented in Fig. 1111 In agreement with 
known theory for vortex lattices, the hexatic symmetry 
is a very robust phenomenon. 



V. CLASSICAL MELTING 

The Bosc-Einstein temperature for a two-dimensional 
noninteracting Bose gas in a harmonic trap is given by 

k B T c =u^C,{2)/N, (58) 

where £(x) is the Riemann zeta function and £(2) ~ 1.28. 
When the temperatures are much lower than this temper- 
ature, there is no thermal cloud present in the gas and 
we can easily extend our analysis to this regime. This 
is experimentally relevant, because the zero-temperature 
limit is difficult to reach [9(j ■ However, according to 
the Mermin- Wagner theorem, a two-dimensional crystal 
cannot exist for nonzero temperatures. For an infinite 
hexagonal vortex lattice this can be seen as follows. Let 
u(x, t) be the displacement field of the vortex lattice. We 
define 

UL((l,t) = {q x u x {q,t) + q y u y (q,t))/q, (59) 
MT(q,*) = {qxUyin, t) - q y u x (<i,t))/q, (60) 

as the longitudinal and transverse fluctuations in the dis- 
placement field, respectively. For long wavelengths the 
action for the vortex lattice is then given by 




(61) 

where n is the average atomic density and c\ and Cqq 
are the compressional modulus and the shear modulus 
of the vortex lattice, respectively. We have used that 
the vortex- vortex interaction decays faster than 1/r 2 , 
since the condensate is assumed to be weakly interact- 
ing. From this action we read off that the dispersion 
of the (almost transverse) Tkachenko mode is quadratic, 
i.e., 

w(q) = 2y/c\CQ & q 2 /n. (62) 

The fluctuations in the vortex positions can then be cal- 
culated as 




which is a converging integral if we realize that the in- 
tegration is over the first Brillouin zone. However, for 
nonzero temperature, we have to add the Bose-Einstein 
distribution l/(e /3a, ' fc ' — 1), which for long wavelengths 
can be approximated by ksT /uj{k). The temperature 
fluctuations are therefore given by 




which diverges logarithmically in the infrared. Therefore, 
the usual Lindemann criterion always predicts melting 
for infinite isolated two-dimensional vortex lattices. For 
a finite system there is a natural infrared cut-off of the 
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divergence, but still there are large contributions from 
low-lying modes. Any finite amount of tunneling, how- 
ever, will turn the system into a three-dimensional sys- 
tem, where the Mermin- Wagner theorem does not apply 
any more. If k denotes the momentum in the z direction, 
we can then define 

tti(q, = {q x u x (q,k,t) + q y u y (q,k,t)/q, (65) 
u T (q, M) = (q x u y (q,k,t) - q y u x (q,k,t))/q. (66) 

Tunneling results in an additional term in the action that 
for long wavelengths is given by 



dt 



dk 
2^ 



Jk 2 (\u L \ 2 + \u T f 



where J is an effective hopping parameter. This still gives 
a quadratic dispersion, but with an anisotropic mass, i.e., 



w(q, k) = 2^(cig 2 + Jfc 2 )(c 66 ? 2 + Jfc 2 )/rL 



(67) 



Writing p = (q Xl q y ,k) and using spherical coordinates 
the dispersion becomes 



uj(p) = 2p 2 \J (ci sin 9 + J cos 2 9) (cqq sin 9 + J cos 2 9) / n. 

(68) 

The fluctuations can therefore be calculated as 

(69) 



(u 2 ) T = k B T 



(ci sin 2 9 + J cos 2 9) (c66 sin 2 9 + J cos 2 
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(2tt) 2 2np 2 
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ci sin 2 9 + J cos 2 9 
cqq sin 2 9 + J cos 2 9 



I c 6 6 sin 2 9 + J cos 2 i 
ci sin 2 9 + J cos 2 i 



which clearly converges in the infrared and the fluctu- 
ations remain finite. For a system of coupled pancake 
Bose-Einstein condensates we can therefore still use the 
usual Lindemann criterion. 

In Fig. E| the vortex-crystal radius is plotted as 
a function of the rotation frequency for a fixed and 
strong tunneling and various temperatures. In the regime 
n/uj± € [.938, .941] there is a dynamical instability to- 
wards elliptic shape deformation. This is indicated by 
a shaded region in the plotted figures. We did not cal- 
culate the fluctuations in this regime. This is related to 
the elliptic shape deformation that occurs before a single 
vortex enters the condensate, and that has been inves- 
tigated before theoretically [77], uR uH MM an d has also 
been observed [sij . Because the unstable mode crosses 
zero, this causes huge fluctuations in the neighborhood 
of this instability. In Fig. ^|we calculate the tempera- 
ture for which quantum fluctuations of the vortex crystal 
are equal to the temperature fluctuations, which defines 
the crossover temperature between quantum and classi- 
cal melting. The temperature should be chosen much 




FIG. 12: Crystal radius R CI normalized to the condensate 
radius R as a function of the rotation frequency for N = 250, 
U = 10,t = 1/20, and T = T c /40 (dotted line), T = T c /50 
(dashed line), and T = T c /100 (solid line). In the shaded 
region there is a dynamical instability towards elleptic shape 
deformation. In the inset the frequency of the unstable mode 
is plotted. 
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FIG. 13: Inverse of the crossover temperature T cross a function 
of fl/uj ± for N = 250 and U = 10. Lines are for t = 1/20 
(solid line) and t = (dashed line). As in Fig. 1121 the shaded 
region is excluded because of the presence of a dynamical 
instability. 



lower to observe the effects of quantum melting. In Fig. 
im the vortcx-crystal radius is plotted as a function of 
temperature for a fixed rotation frequency. 

When tunneling between sites is suppressed, the cor- 
relation between neighboring vortices 



((Ui„ - u„„) 2 ) 



2A 2 



(70) 



where A mn is the distance between the vortex n and 
m at site i, has been proposed as an appropriate order 
parameter, with unchanged Lindemann parameter [49l 
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FIG. 14: Crystal radius R CI normalized to the condensate 
radius R as a function of T c /T for N = 250, U = 10, and 
Q/co± = 0.93. From left to right the lines have hopping pa- 
rameter t = 1, t = 0.1, t = 0.05, and t = 0.01, respectively. 

l9ll 193 . 193^ . In this order parameter the diverging small 
momentum contributions are subtracted. Nevertheless, 
it turns out that we have to go to low temperatures to 
see any crystalline order. The same phases as in the case 
of zero temperature can be distinguished. The result of 
this calculation is plotted in Fig. ^] Note that a true 
solid phase of the vortex lattice is not observable in these 
phase diagrams, not even for the lowest temperature of 
T = T c /100. 



VI. ANHARMONIC RADIAL CONFINEMENT 

The inclusion of a quartic potential, in addition to the 
usual harmonic potential, has attracted interest for sev- 
eral reasons. The quartic potential plays an important 
role in the understanding of vortex nucleation HH • 
For our purposes, however, it is more important that in 
this way rotation frequencies larger than the trap fre- 
quency can be applied. This causes the density to be 
lower in the center of the trap, and eventually gives rise 
to the formation of giant vortices, i.e., multiple quantized 
vortices, in the center of the trap j5j| |5(| [§3 Experi- 
ments in this setup were performed and it was observed 
that the vortex lattice became disordered, but no giant 
vortex was observed [l5j . Giant vortex formation was 
observed by artificially removing the inner part of a fast 
rotating condensate by means of a laser beam [TlT | . Shape 
oscillations of a vortex lattice in an anharmonic poten- 
tial were also studied experimentally j9||. Theoretically, 
the phase diagram of this system was studied intensively 
to identify the parameter space where giant vortex for- 
mation can take place IH, • Other aspects that 
are studie d are the dynamics of forming th e gia nt vortex 
jlOOt llOlj , oscil lations of the vortex lattice |l02j| , aspects 
on observation |l03t Il04j , and stability of quantum fluc- 



13 




0.92 0.93 0.94 0.95 0.96 0.97 0.98 




0.92 0.93 0.94 0.95 0.96 0.97 0.98 




0.92 0.93 0.94 0.95 0.96 0.97 0.98 



FIG. 15: Phase boundaries between the solid of solid shells 
(SS), hexatic liquid (HL) and liquid (L), as a function of radial 
distance and rotation frequency for N = 250, U = 10, and 
t = (dashed line) and t = 1 /20 (solid lines) . The figures are 
for a) T = T c /20, b) T = T c /35, and c) T = T c /100. 

tuations |l05l | . 

It is straightforward to extend our analysis to this case. 
We first of all add the quartic potential 

V i {r)=\^, (71) 
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FIG. 16: (Color online) Classical vortex lattice and density 
profile for rotation frequencies Sl/cu± = 1.05, 1.1, 1.5 in 
the presence of a quartic potential. Parameters are chosen as 
U — 10 and A = 0.01. White means high density, black low 
density. The vortex positions are indicated by a dot, such 
that also the vortices outside the condensate are visible. 



FIG. 17: Inner and outer radius of the condensate as a func- 
tion of the rotation frequency for N = 250, U = 10, A = 0.01, 
and t — 0. The solid line is the result of a calculation using the 
lowest Landau approximation and with 37 vortices taken into 
account. The dashed line is the result of the approximation 
that the density profile is described by a quartic polynomial. 

where £± is the harmonic length associated with the ra- 
dial trapping. We again use lowest Landau level wave- 
functions and make the ansatz that there is one vortex in 
the center and that the rest of the vortices order them- 
selves in rings of multiples of six. The resulting density 
profiles for 37 vortices are plotted in Fig. ^|as a func- 
tion of the rotation frequency. We compare these density 
profiles with the Thomas-Fcrmi-likc density profiles of 
the form 

^° c (5 ±i )(4" 1 )' (72) 

that come from minimizing the on-site energy functional 
without paying attention to the lowest Landau level con- 
straint. When the minus sign is chosen there appears a 
hole in the density, and the condensate shape is annular 
with inner radius R\ and outer radius i?2- We compare 
the values coming from this ansatz with the inner and 
outer radius coming from the lowest Landau level densi- 
ties in Fig. El They agree quite well. 

We also calculate the quantum fluctuations of the vor- 
tices. This is only possible in a limited regime, where 
the ansatz of a central vortex surrounded by rings of six 
vortices is dynamically stable. In this regime, the inner 
part of the vortex crystal is melted, although the particle 
density is nonzero there. We define the liquid radius i?L 
as the radius of the innermost vortex that is part of the 
vortex crystal. In the region fl/uj± € [1.155,1.25] the 
ansatz is stable against fluctuations. We find that the 
liquid radius in this region is almost constant and given 
by Rl/R — -58. There is a second radius which separates 
the vortex crystal from the vortex liquid that is at the 
outside of the condensate, which corresponds to the crys- 
tal radius R cl - defined before. In this region this radius 
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is always equal to the condensate radius or R cr /R = 1.0. 

VII. CONCLUSIONS AND OUTLOOK 

In this article we derived the theory of vortex fluctu- 
ations in a one-dimensional optical lattice. We showed 
that in this configuration the modes of the vortex lat- 
tice get a dispersion in the axial direction. In particular 
we discussed the Bloch bands of the Tkachenko modes. 
Based on the Lindcmann criterion we studied the melting 
of the vortex lattice. Because of the inhomogeneous den- 
sity the melting occurs from the outside inwards. Com- 
parison with a local density-approximation yields impor- 
tant finite-size corrections, because the local-density the- 
ory does not take into account the discrete nature of the 
vortices and overestimates the rotation frequency that is 
needed for total melting of the vortex lattice. Tunneling 
between the sites of the optical lattices decreases the fluc- 
tuations and causes freezing of the vortex lattice. Tem- 
peratures on the other hand increases the fluctuations 
considerably. The crossover temperature from classical 
to thermal melting is very low. which makes it an exper- 
imentally difficult problem to see the effects of quantum 
fluctuations. By looking into the correlations between 
the vortices, several phases can be distinguished, where 
part of the order is destroyed by quantum fluctuations. 

In the one-dimensional optical lattice there is a clear 
experimental signal for both the fluctuations in the vor- 
tex position and the liquid. The fluctuations can be mea- 
sured in analogy with the situation of a single vortex in 



an optical lattice. By imaging in axial direction one will 
see the vortex cores distributed themselves in a gaus- 
sian distribution around their equilibrium position, from 
which the size of the fluctuations can be extracted [4(| ■ In 
the liquid the vortices are no longer individually visible, 
which is a clear distinction from the vortex crystal. An 
interesting question is whether the liquid will completely 
restore the rotation symmetry that is broken by the pres- 
ence of the vortex lattice, or the liquid is partly pinned 
because of the interaction between the vortex liquid and 
the vortex crystal. 

It remains a challenging problem to describe the co- 
existing crystal-liquid. This will allow to decide on the 
occurrence of melting based on energy considerations and 
thus shed more light on the accuracy of the application 
of the Lindcmann criterion in this inhomogeneous situa- 
tion. This also applies to the phases where partial melt- 
ing takes place. It is not only important to find a good 
description of these phases themselves, but also the pre- 
dicted phase coexistence between them raises interesting 
new questions. 
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